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Three interaction site (TIS) model of RNA 

We develop a realistic force field for nucleic acids using the TIS modelPin which each nucleotide 
is replaced by three spherical beads P, S and B, representing respectively a phosphate, a sugar and 
a base (Figure SI). The coarse-grained beads are at the center of mass of the chemical groups and 
have a radius Ri of 2 A for phosphates, 2.9 A for sugars, 2.8 A for adenines, 3 A for guanines 
and 2.7 A for cytosines and uracils. The values of Ri are calculated using Vi = AttR] /3, where 
is the van der Waals volume of the chemical group computed from the coordinates and radii of 
its individual atoms. We use the total molecular weight of each RNA group as the mass of the 
representative bead in our simulations. In the TIS representation of nucleic acids, bond lengths, p, 
and valence angles, a, are constrained by harmonic potentials, U (p) = kp{p — po)^ and U (a) = 
ka{(X, — Oof- , where the equilibrium values po and ccq are obtained by coarse-graining an ideal 
A-form RNA helix.'^The values of kp, in kcalmol^^A^^, are: 64 for an S(5')— P bond, 23 for an 
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P— S(3') bond and 10 for an S— B bond. The values of ka are 5 kcalmol^^rad^^ if the valence angle 
involves a base, and 20 kcalmol^rad^ otherwise. We have chosen kp and ka so that the time 
averages of (p — po)^ and (a — ao)^ measured in simulations at 15 °C, match the corresponding 
quantities averaged over all bonds in the coarse-grained NMR structure of the hTR pseudoknot 
(PDB code 2K96). 

Single strand stacking interactions, U^i, are applied to all pairs of consecutive nucleotides 
along the chain, 

^ST 

" 1 + 1.4(r - ro)2 + 4(01 - 0io)2 + 4(02 - 02o)2 ' 
where r, 0i and 02 are defined in Figure SI a. The equilibrium values tq, 0io and 020 are extracted 
from the coarse-grained structure of an ideal A-form RNA helix ^ and depend on the chemical 
identities of the two nucleotides. We obtain the constants U^r^ for sixteen distinct nucleotide dimers 
from available experimental data on stacking of nucleic acid bases in single- stranded and double- 
stranded RNA,'^^as described next. 

Additive contributions AG(^I^) of individual stacks to the total stability of an RNA double 
helix, where a — d and b — c are stacked Watson-Crick base pairs, are known experimentally.'^ For 
reference, experimentally determined enthalpic AH{^Z^) and entropic A5(*l'y) contributions to 
AG(*I^) are reproduced in Table SI. We make the following approximations: 

^(^a^d) " AH(^^+AH(^^+0.5AH{a-d)+0.5AH{b-c), 
AsC-') = AsO+Asf"), (2) 



a — dj \aj \c 

where AH (^) and A5(^) are the enthalpy and entropy changes resulting from stacking of b over a 
along 5' — 7- 3' in one strand, respectively, and AH (a — d) is the additional stability due to hydrogen 
bonding between a and d in two complementary strands. Inspection of A// (^3^) and A5(^I^) in 
Table SI leads us to conclude that, except for (q) and (^), stacking parameters do not depend 
strongly on the order of nucleotides along 5' — )• 3'. Therefore, we assume, with the exception of 
Q and (g, that AhQ = AH{1) and A5(*) = A5(^), which is also valid within the range of 
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experimental uncertainties specified in ref 2. This assumption allows us to average the experimen- 
tal values for A/zgl^) and AH(£:^), AH{^Zg) and AH{^Z^), AH{^Zg) and AH{^Z^), and 
similarly for the corresponding entropies. 

Based on the experimental data for stacking of nucleic acids in single- stranded RNA (see Table 
8. 1 in ref 3 and Table 1 in ref 4), we make additional assumptions that AH (^) = AH (^) , AH (^) = 
A//(^) = A//(^), and similarly for the entropies. This allows us to combine eq 2 for (^ly) with 
the experimentally determined melting temperature (^) = 26 °C,^ which under our assumptions 
equals Q , and to solve for AH Q , AS Q and A//( A - U) . Now putting AH{^)=AH Q and 
A5(^) = A5(^) in eq 2 for (^1^) , we can compute AH (y) and A5(y) . Finally, we assume 

As(^^j =tAs(^")+(l-t)As([^j, (3) 

where ^ = 0.615 yields ?m(^) = 13 °C, which matches the experimental result for t,n{^^Thc 
remaining stacking parameters follow directly from eq 2 without any additional approximations, 
if we use the computed hydrogen bond enthalpy AH (A — U) = — 1 .47 kcal/mol for an A — U base 
pair and 3/2 times this value for a G — C pair. The resulting stacking parameters, which are used 
in the present simulations, are given in Table S2. The relative stabilities of stacks agree with the 
experimental data,'^ identifying tm (q) and tm (jj) as the highest and lowest melting temperatures 
among all stacks. 

We simulated stacking of nucleotide dimers, similar to that shown in Figure SI a, using the 
stacking potential Usj in eq 1 and U^j = —h -\-k^{T — Tm)s, where is the Boltzmann constant, 
T (K) is the absolute temperature, Tm (K) is the melting temperature of each stack (from Table S2) 
and h and s are adjustable parameters. In our simulations, we computed the stability AG{T) of 
stacks at temperature T as 

AG{T) = -kBT\ogNi+kBTlogN2 + AGo, (4) 
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where A'^i is the number of all stacked configurations for which Ust < —k^T and N2 is the number 
of all unstacked configurations. We adjusted h and s individually for all stacks so that the simu- 
lation results for AH and AS, given by AG{T) = AH — TAS, matched the corresponding values in 
Table S2. The stability correction AGq in eq 4 is assumed to be constant for all stacks and accounts 
for potential discrepancies between AG{T) measured in experiments and its definition used in our 
simulations. 

As an example, in Figure S2 we plot AG(T) obtained in simulations of the (^) stack with 
various values of h and s = 0. The simulation results in Figure S2 are shown for AGq = 0. The 
observed melting temperature T*, defined as AG{T*) = 0, increases with h and equals the target 
melting temperature 7m (^) in Table S2 when h = 5.98 kcal/mol. For s = the entropy loss A5 
associated with stack formation, given by the slope of AG{T) over T, is smaller than the value of 
A5(^) specified in Table S2. To correct this, we use U^^ = —5.9S + kB{T — Tm)s, which does not 
result in changes in the melting temperature but allows us to adjust the slope of AG{T) by adjusting 
the value of s. We find that s = 5.30 yields AS{^) in Table S2. 

We carried out the same fitting procedure for all nucleotide dimers. The resulting parameters 
U^j for AGo = and AGq = 0.5 kcal/mol are summarized in Table S3. In our simulations, we use 
the value AGq = 0.5 kcal/mol which yields the best agreement with experiments (see discussion in 
the last section below). Note that, although some stacks have equivalent thermodynamic parame- 
ters in Table S2, they may require somewhat different U^j due to their geometrical differences. 

Coarse-grained hydrogen bond interactions t/nB are assigned based on the hydrogen bonds 
present in the original NMR structure. In this work, we carried out independent simulations of 
the pseudoknot and hairpin conformations of the hTR pseudoknot domain (PDB codes 2K96 and 
1NA2, respectively). In both cases, we generated an optimal network of hydrogen bonds by sub- 
mitting the NMR structure to the WHAT IF server at| |http : //swift . cmbi . ru . nl| Each of 
the generated bonds is modeled by a coarse-grained interaction potential, 

Uhb = 2.286 X [1 + 5(r - rof + 1.5(0i - Oio? + 1-5(02 - 02o)^ 



S4 



+0.15(v/-V/o)^ + 0.15(v/i-V/io)2 + 0.15(v/2-V^2o)^] \ (5) 



where r, Oi, O2, V^, V^i and i//2 for various coarse-grained sites are defined in Figure SI. In the case 
of Watson-Crick base pairs, the equilibrium values ro, 0io, 020^ V^o^ V^io and 1/^20 are adopted from 
the coarse-grained structure of an ideal A-form RNA helix.l^For all other bonds, the equilibrium 
parameters are obtained by coarse-graining the PDB structure itself. Equation 5 specifies Uub for a 
single hydrogen bond and it must be multiplied by a factor of 2 or 3 if the same coarse-grained sites 
are connected by more than one hydrogen bond. The complex geometry of Uub is the minimum 
necessary to maintain stable double (and triple) helices in our coarse-grained model. 

Comparison of two alternative sets of interaction parameters 

In our simulations, stacking parameters U^j are determined based on the definition of AG{T) given 
in eq 4. The corrective constant AGq = 0.5 kcal/mol in eq 4 is introduced to improve quantitative 
agreement between simulation and experimental melting data of the hTR pseudoknot domain. If 
AGq were omitted from eq 4, this would result in stronger stacking interactions (see Table S3). 
In order to preserve overall melting temperatures, this increase in the magnitude of t/g-p must be 
compensated for by a decrease in the strength of hydrogen bonds. For AGq = 0, the best agreement 
with experiments is achieved when the prefactor in eq 5 is reduced to 2.065. 

In Figure S3 we compare melting data for the hairpin (HP) conformation of the hTR pseudo- 
knot domain for the two parameter sets, with and without the corrective constant AGq. The melting 
profile of the Watson-Crick part of the double helix, stem 1 in HP, is hardly affected by the choice 
of U^j. However, we observe a large discrepancy for the uridine-rich stem 2, whose melting tem- 
perature increases from 50 °C to 75 °C if AGq is set to in eq 4. In the latter case, the total melting 
profile of HP shows only one peak at 78 °C, which does not compare well with two experimental 
peaks at 50 °C and 79 °C.^ We conclude that the corrective constant in eq 4 is crucial for obtaining 
quantitative agreement with the experimental data for HP. 



S5 



Melting of secondary structure in the pseudoknot (PK) conformation of the hTR pseudoknot 
domain is compared in Figure S4 for AGq = 0.5 kcal/mol and AGq = 0. In experiments,'^ the 
temperature range for melting of stems 1 and 2 in PK is 65-95 °C and 60-80 °C, respectively. As 
shown in Figure S4, these experimental data are reproduced well in simulations with AGq = 0.5 
kcal/mol (weak stacks in Table S3 and U^f^ = 2.286). At the same time, for AGq = (strong stacks 
in Table S3 and = 2.065), melting of both stems occurs in the temperature range 50-95 °C. 
For both parameter sets, the overall melting profile of PK has a sharp peak at 70 °C, in agreement 
with experiments.'^ 

The distance between two peaks in the melting profile of HP increases with AGq and exceeds 
the experimental distance when AGq > 0.5 kcal/mol (data not shown). We therefore conclude that, 
for both conformations of the hTR pseudoknot domain, the parameter model based on AGq = 0.5 
kcal/mol yields the optimal agreement with experimental thermodynamic data. 



Crowder-RNA interactions 



Crowder-RNA interactions are modeled by a generalized Lennard- Jones potential. 



Uu{r) 
Uu{r) 



2Ri 



r + Do-D 



12 



r + Do-D 



+ 1 



0, r>D, 



(6) 



where r is the distance between the particles' centers of mass, Dq = 3.2 A is the effective penetra- 
tion depth, Ri is the radius of an RNA coarse-grained bead (values specified above), rc is the radius 
of a crowder, and D = Ri + rc- The ratio IRj/DQ in eq 6 is used to scale the interaction strength 
e = 1 kcal/mol in proportion to the surface contact area. 

We use the same formula to model RNA-RNA excluded volume interactions, but take Ri = 1.6 
A for all RNA beads. In this case, eq 6 becomes a standard (purely repulsive) Lennard- Jones 
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potential. 



Uu(r) = .[(^) 



12 




6 



-2 



+ 1 , r < Do, 



f/uW = 0, r>Do. 



(7) 



With adjustment of Rj steric claslies between two stacked bases are avoided. 

Simulation details 

The RNA and crowder dynamics are simulated by solving the Langevin equation, which for par- 
ticle / is rriifi = —YiTi + Vi + R(, where m/ is the particle mass, Ji is the drag coefficient, F/ is the 
conservative force, and R, is the Gaussian random force, (R,(?)Rj(?')) = 6k^Tyi5ij5{t — t'). The 
drag coefficient is given by the Stokes formula, Ji = 6Tzr]Ri, where r\ is the viscosity of the medium 
and Ri is the particle radius. To enhance conformational sampling, we take r\ = 10~^Pa s, which 
equals approximately 1% of the viscosity of water. The masses and radii of RNA coarse-grained 
beads are specified above. The masses of crowders scale with their volume, assuming density equal 
to that of a typical folded protein such as ubiquitin, which has the molecular weight Mr = 8564 Da 
and radius rc = 1.2 nm. The Langevin equation is integrated using the leap-frog algorithm with 
time step A? = 2.5 fs. The length of a simulation run at each temperature is 2.5 /is. 

The number of crowders of a given type is computed from its specified volume fraction and 
the volume of the simulation box. In simulations with large crowders (rc = 10.4 nm, 5.2 nm and 
2.6 nm) we use 60 nm as the side of the cubic simulation box. For example, the E. coli mixture 
contains the volume fractions 0i = 0.11, 02 = 0.11 and 03 = 0.08 of crowders with rc = 10.4 nm, 
5.2 nm and 2.6 nm, respectively. For a simulation box with side 60 nm, this yields 5 crowders 
with rc = 10.4 nm, 40 with rc = 5.2 nm and 234 with rc = 2.6 nm (279 crowders in total). In 
simulations with small crowders (rc =1.2 nm and 0.6 nm) the number of crowders in the cubic 
box with side 60 nm becomes very large. To minimize simulation time of these systems, the 
size and shape of the simulation box is adjusted periodically to accommodate RNA in its current 
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conformation. To do so, we first align the walls of the box with the RNA axes of inertia. The new 
position of the RNA center of mass and new side lengths L^, Ly, are computed so that there is 
at least a 6 nm distance between each RNA bead and all six walls. When, as a result of diffusion, 
the distance between an RNA bead and a wall becomes less than 2.4 nm, the box is adjusted again. 
The size and shape of the simulation box are therefore directly coupled to the RNA configuration. 
For instance, in simulations at °C, the RNA remains folded and L^, Ly, fluctuate around 18 
nm, 15 nm and 14 nm, respectively (assuming > Ly > L^). In simulations at 120 °C, the RNA 
is unfolded and the average lengths are L^ = 20 nm, Ly = 16.5 nm and = 15 nm. In simulations 
at intermediate temperatures, when the RNA folds and unfolds, the size of the box undergoes large 
fluctutions between the high-temperature and low-temperature values. The frequency with which 
the simulation box is adjusted also depends on the temperature through RNA diffusion. The box 
is adjusted approximately every 1200000 steps at °C and every 700000 steps at 120 °C. When 
the box is adjusted, the number of crowders changes with the new box volume in order to keep the 
volume fractions ^ constant. 
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Table S 1 : Theraiodynamic parameters of double-stranded stacks from ref 2. In the first column, 
the 5' to 3' direction is shown by an arrow. 





AH, kcalmol""^ 


AS, calmol~^K~^ 


A— U 
A-U 


-6.82 


-19.0 


U— A 
A-U 


-9.38 


-26.7 


A— U 
U-A 


-7.69 


-20.5 


U-A 
C-G 


-10.48 


-27.1 


A— U 
C-G 


-10.44 


-26.9 


T T A 

U — A- 

G-C 


-11.40 


-29.5 


A-U 
G C 


-12.44 


-32.5 


G-C 
C-G 


-10.64 


-26.7 


G-C 
G-C 


-13.39 


-32.7 


C G 
G-C 


-14.88 


-36.9 
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Table S2: Thermodynamic parameters of single- stranded stacks, derived in this work. The match- 
ing enthalpies of hydrogen bond formation in Watson-Crick base pairs are given in last two rows. 
The melting temperatures of stacks are indicated in °C, tm{°C) = Tm(K) — 273.15. In the first 
column, the 5' to 3' direction is shown by an arrow. 



1 a 


AH, kcalmol ^ 


AS, calmol-iK-i 


t °C 


U 

u 


-1.81 


-7.2 


-21 


c 
c 


-2.87 


-10.0 


13 


c. u 
U' c 


-2.87 


-10.0 


13 


A 
A 


-3.53 


-11.8 


26 


A. U 
U' A 


-3.53 


-11.8 


26 


A. C 

C A 


-3.53 


-11.8 


26 


G 
C 


-4.21 


-13.3 


42 


G. U 
U' G 


-5.55 


-16.4 


65 


C 
G 


-6.33 


-18.4 


70 


G. A 
A' G 


-6.75 


-19.8 


68 


G 
G 


-8.31 


-22.7 


93 



AH {A - U) = - 1 .47 kcal/mol 
AH{G-C) = -2.21 kcaymol 
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Table S3: Temperature-dependent stacking parameters U^j used in eq 1. The two sets of h corre- 
spond to two different values of the additive constant AGq in eq 4: 0.5 kcal/mol and (numbers in 
brackets). The melting temperatures of individual stacks are given in Table S2. The values of s 
conform to kB{T — Tm) evaluated in kcalmol"^ In the first column, the 5' to 3' direction is shown 
by an arrow. 





Ul^=-h + kB{T-Trr,)s 


' a 


h, kcalmol ^ 


s 


U 
U 


3.52 (4.27) 


-3.56 


C 

c 


4.16 (4.87) 


-1.57 


c. u 
U' c 


4.14 (4.88); 4.14 (4.87) 


-1.57; -1.57 


A 
A 


4.49 (5.19) 


-0.32 


A. U 
U' A 


4.43 (5.16); 4.45 (5.15) 


-0.32; -0.32 


A. C 

C A 


4.43 (5.16); 4.45 (5.15) 


-0.32; -0.32 


G 
C 


4.75 (5.48) 


0.77 


G. U 

U' G 


5.17 (5.89); 5.12 (5.84) 


2.92; 2.92 


C 
G 


5.22 (5.93) 


4.37 


G. A 
A' G 


5.26 (5.98); 5.22 (5.95) 


5.30; 5.30 


G 
G 


5.70 (6.42) 


7.35 
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(a) (b) 




Figure SI: Geometrical parameters for (a) stacking in eq 1 and (c-d) hydrogen bonding in eq 5. 
Sites P, S and B are shown in black, green and red, respectively. Sample conformations are derived 
from 2K96.pdb. The numbers refer to specific nucleotides, and r, 0, and or i// refer to distances 
(dist.), bond angles (ang.) and dihedral angles (dih.) between indicated sites, (a) Stacking: r = 
dist. (B102, B103), (^i = dih. (P102, S102, P103, S103), (^2 = dih. (P104, S103, P103, S102). (b) 
Hydrogen bonding between B and B: r = dist. (B98, B172), 0i = ang. (S172, B172, B98), 62 = 
ang. (S98,B98,B172), v/ = dih. (S98, B98, B172, S172), i/Zi = dih. (B98, B172, S172, P173), 1^2 
= dih. (B172, B98, S98, P99). (c) Hydrogen bonding between B and S: r = dist. (S98, B172), 0i 
= ang. (S172, B172, S98), 02 = ang. (P99, S98, B172), i// = dih. (P99, S98, B172, S172), i/Zi = 
dih. (S98, B172, S172, P173), V/2 = dih. (B172, S98, P99, S99). (d) Hydrogen bonding between P 
and S: r = dist. (S102, P176), 0i = ang. (S176, P176, S102), 02 = ang. (P103, S102, P176), i/a = 
dih. (P103, S102, P176, S176), i^i = dih. (S102, P176, S176, P178), ^2 = dih. (P176, S102, P103, 
S103). 



S13 



0.60 0.64 0.68 0.72 0.76 
k^T, kcal / mol 

Figure S2: Stability AG(r) of stack (^) , computed from eq 4 with AGq = 0, in simulations with the 
stacking potential given by eq 1 and f/g^ = —h + k^{T — Tm)s. indicates the melting temperature 
for (^) from Table S2, k^T^ = 0.68 kcal/mol. Open and closed symbols show simulation results 
for various h and s. Red solid line indicates the target stability AG{T) = AH — TAS for (^) , where 
AH and AS are given in Table S2. Same stability line is obtained in simulation with h = 5.98 
kcal/mol and = 5.30 (closed symbols). 
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Figure S3: (a) Number of intact base pairs A^bp vs. temperature in two elements of secondary 
structure in HP. (b) The melting profiles show the rate of change of the number of intact base pairs 
with temperature. The data are obtained for different values of AGq in eq 4: 0.5 kcal/mol (solid) 
and (dashed). 
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Figure S4: Same as in Figure S3, but for the PK conformation. 
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Figure S5: Temperature dependence of the stabilities of the hairpin and pseudoknot in the absence 
of crowders (solid) and in a monodisperse suspension of crowders with ^ = 0.3 and rc = 1.2 nm 
(dashed), (a) AU177 sequence, (b) AU177 sequence with DKC mutations. 
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Figure S6: (a) Changes in stability (kcal/mol) of HP and PK at 37 °C due to crowders, as a function 
of the crowder radius re- (b) Probability distributions p{Rg) of the radius of gyration of the 
unfolded HP and PK structures at 37 °C in the absence of crowders (0 =0). Rq of HP is defined 
as the radius of gyration of strand G93-C121, which is the only fragment that is structured in HP 
(same definition is used in the paper). Rq of PK is computed for the entire length of RNA. The 
unfolded state cannot be sampled directly in simulations at 37 °C and the shown p{Rg) are obtained 
by statistical reweighting of the high-temperature data. For PK, p{Rg) is noticeably asymmetric, 
showing an extended tail on the side of large Rq. This is in part due to residual stacking in the 
unfolded state, whose effects become more pronounced for longer strands. A larger mean of p{Rg) 
and its asymmetry explain why the PK structure is more strongly influenced by crowders than HP. 
(c)-(d) p{Rg) in crowder suspensions with ^ = 0.3 and different re- Symbols match those in 
panel (a) and black solid curves show p{Rg) at ^ = 0. In all four panels, dashed curves are drawn 
through data to facilitate comparison. 
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Figure S7: Same as in Figure 2a of the main text, but in a monodisperse suspension of crowders 
with = 0.3 and rQ — l.l nm. 
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Figure S8: Melting profiles of the hairpin in the absence of crowders (black) and in a monodisperse 
suspension of crowders with = 0.3 and rc = 1.2 nm. Dashed and solid curves are for the AU177 
sequence with and without additional DKC mutations. 
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